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Abstract 

We describe how to evaluate approximately various physical interesting quantities in random 
Ising systems by direct renormalization of a finite system. The renormalization procedure is 
used to reduce the number of degrees of freedom to a number that is small enough, enabling 
direct summing over the surviving spins. This procedure can be used to obtain averages of 
functions of the surviving spins. We show how to evaluate averages that involve spins that 
do not survive the renormalization procedure. We show, for the random field Ising model, 



how to obtain T{r) = {a{0)a{ r )) — {a{0)){a{ r )) , the "connected" correlation function and 



S{r ) = (o"(0)(t( r )), the "disconnected" correlation function. Consequently, we show how to 
obtain the average susceptibility and the average energy. For an Ising system with random bonds 
and random fields we show how to obtain the average specific heat. We conclude by presenting 



our numerical results for the average susceptibility and the function r(r) along one of the principal 
axes. (We believe this to be the first time, where the full three dimensional correlation is calculated 
and not just parameters like or r].) The results for the average susceptibility are used to extract 
the critical temperature and critical exponents of the 3D random field Ising system. 

PACS numbers: 05.50.+q, 64.60.Cn, 75.10.Nr, 75.10.Hk 
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I. INTRODUCTION AND OUTLINE 



Real space renormalization group (RSRG) served as a major tool, over the last thirty 
years, in the field of critical phenomena. By simplifying calculations near a critical point, 
the various RSRG techniques, such as using the majority rule Q|, the well known Migdal- 
Kadanoff (MK) , the Casher-Schwartz and others [5|, allow one to penetrate the 

critical regime to a point at which critical exponents can be extracted. For translational in- 
variant (pure) systems, only a single renormalization step is required to obtain the recursion 
relations for the parameters of the Hamiltonian, from which critical fixed-points may be de- 
rived along with critical exponents. For random systems the recursion relations are position 
dependent. Therefore, a natural approach is to consider the recursion for the distribution of 
disorder. Equivalently, recursion relations may be obtained for all the parameters defining 
the distribution, moments, correlations, etc. Practically, in this approach, the recursion 
relations are truncated to obtain relations involving only the mean and variance, keeping 
the random couplings^dependent BBBB. An alternafve approach, suggested first by 
Berker and Ostlund [10||, is to consider a given realization of disorder on a finite system. 
Renormalization is then used to reduce the system to a size where brute force calculation is 
possible. Thermal averages of certain quantities can thus be obtained for that realization and 
ensemble average is obtained by repeating the procedure for many realizations and averaging 
the results. The advantage of the method is that all the moments and correlations generated 
by renormalization are kept. The disadvantage is that the renormalization leaves in the end 
a small number of spins and, therefore, only thermal averages of functions of those spins can 
be evaluated directly. This is good enough to obtain directly the ensemble average of the 
magnetization because the average magnetization obtained from the surviving spins is 
exactly the true average magnetization. If, on the other hand, we are interested in averages 
involving spins that do not survive the renormalization process, things become much more 



complicated. Take for example the ensemble averaged correlation, T{rij) = {(Ticrj) — {ai){aj), 
where (■ ■ ■ ) denotes thermal average and denotes ensemble average. It can be calculated 
directly from the remaining spins provided that the vector rjj connecting the sites i and j 
equals a vector connecting two of the surviving spins or obtained from it by a symmetry 
operation on the initial lattice. In any other case, a direct calculation is impossible. It is true 
that quantities like the susceptibility that involves Fj^'s for all pairs of sites can be calculated 
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indirectly 
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1 Following Dayan et al jlii|], this requires adding a constant field, 

H , to the system, calculating the ensemble average of the magnetization as a function of H 
by the method outlined above and then differentiating with respect io H ai H = In fact, 
many interesting quantities may be obtained by adding the appropriate interaction to the 
Hamiltonian and differentiating ensemble averages with respect to the corresponding cou- 
pling constants. The trouble with that approach is that even in the relatively simple case of 
evaluating the susceptibility, the numerical differentiating is quite problematic. Berker and 
coworkers were using the chain rule to approximately recover thermodynamic 

densities of the original system from the renormalized couplings of the reduced system. The 
main problem here is that the method is limited only to the obtainment of thermal averages 
of products of spins showing in the Hamiltonian. The purpose of the present article is to 
show how to calculate various interesting quantities that involve spins that do not survive 
the renormalization by a direct and effective method. 

The paper is organized as follows. In Sec. m we describe briefly the Casher Schwartz 
(CS) renormalization procedure, which is the RSRG that we use here for the numerical 
demonstration of our method in Sec. IIVI Note, however, that the method we present is 
more general and can be used with any other renormalization scheme. In Sec. IIIIl the 
elements of our method are mainly considered for the random field Ising system. It is shown 
how to calculate the average "connected" spin-spin correlation function, and the average 



"disconnected" spin-spin correlations, Sij = {(JiCTj), from which the average susceptibility, X: 
and the average total energy, E, can be easily obtained. It is further shown how to calculate 
the average specific heat, C, for an Ising system with random bonds and random fields. Note 
that the method presented here enables a full evaluation (though approximate) of quantities 
like Tij and Sij that depend on distance. As we shall show later, this follows from two facts: 
(a) The finite and small number of spins we are left with at the end of the renormalization 
procedure, (b) The fact that the system is random. (It will be shown how, in principle, this 
method can be used to calculate by a similar method of a finite system renormalization, 
such quantities in the pure system. The practicality of the method for the pure system will 
prove, however, to be questionable.) In the last section, we demonstrate the usefulness of 



our method by calculating the average susceptibility x and F(r) for r's lying on a main axis 
of the lattice, for the random field Ising system. The evaluation of x is used to derive critical 
exponents that may be compared with the exponents derived by other methods. Note that 



3 



the values of the exponents depend not only on the numerical application of the method 
presented here, but also on the specific scheme of renormalization employed. 



II. RENORMALIZATION 

Although our method is general, we will use the CS scheme ^ for the numerical demon- 
stration of our method and present the results in Sec. IIVI Like any other renormalization 
procedure (such as Migdal-Kadanoff (MK) and others jO] ) , when performed on a regular 
lattice, recovering the original form of the Hamiltonian is not an exact procedure. Never- 
theless, for the translational invariant (pure) Ising system it produces good results ^. In 
most RSRG calculations for random systems, correlations generated by the renormalization 
are simply ignored. It was suggested, however, many years ago by Harris and Lubensky Q| 
that those correlations are important. The CS scheme generates, indeed, such correlations. 
Schwartz and Fishman j^, used the CS scheme to renormalize the mean and variance of 
the distribution of random bonds. They took into account the generated correlations and 
found that inclusion of the effect of generated correlations in the the renormalized variance 
is essential. 

The CS renormalization method for the pure Ising system is described in detail in P|, 
while a detailed demonstration of how it can be used, locally, for a random bond system can 

n 

be found in p]. Moreover, our numerical renormalization procedure here, as conducted for 
the random field Ising model, follows almost exactly the numerical procedure used by Dayan 
at al Here, therefore, we only describe it in brief. According to the CS scheme, an 

integration of every other site is performed exactly, but then, all non-nn bonds, generated 
by the RG transformation, are symmetrically bent onto available nn bonds, many-spin odd 
interactions may be grouped to form the renormalized field, while many-spin even inter- 
actions are simply omitted. Here, though, as in Ref. in order to simplify computer 
programming, we only keep the renormalized fields, and ignore the many-spin odd interac- 
tions as well. The integration of every other site is, relatively, an easy task, even in 3D, 
since every spin situated on an even site interacts only with neighboring spins situated on 
odd sites. 

We start then with a set of = Ising spins, cTj = ±1, with L = 2^, situated on 
a 3D SC lattice. Suppose, now, that the Ising system is not translational invariant and 
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represented by the Hamiltonian 

H = - ^ JijOiOj - '^hiCXi, (1) 

<i,j> i 

where < i,j > refers to nearest neighbors only. Performing the trace over every other site, 
each of the erased spins, contributes separately to each of the terms in the new Hamiltonian, 
generating all possible interactions among its 6 nn's. We arrive, then, at a new Hamiltonian, 
containing fields, nn, nnn and multi-spin interaction terms, from 3-spin and up to 6-spin 
interactions. All the couplings are again local. The result is even further complicated by 
the fact that the resulting lattice is not an SC but an FCC lattice. As was mentioned 
above, the multi-spin interactions are simply ignored, while the values of the 3 generated 
nnn interactions are symmetrically distributed over the 12 nn interactions. To bring the 
lattice back to its SC form, we still need to integrate over each of the face-centered spins. To 
do that, we first bend all nn bonds connecting between face-centers, onto nn bonds that lie 
on the face of the cube and connecting between face-centers and vertices. The extra erasure 
step can now be easily executed following the CS 2D renormalization scheme. This is a 
much simpler procedure, which we shell not describe here, but can also be found in Refs. 

m. 

In our study, the above two steps procedure is performed locally and repeated iteratively 
until the system is brought down to a size of 2 x 2 x 2, for which a trace can be performed 
exactly. 

III. THE METHOD 

In this section we describe, mainly, how to calculate the "connected" and "disconnected" 
correlation function for the three dimensional random field system. We consider the random 
field Hamiltonian 

U = - J ^ (Jiffj - ^ hiai. (2) 

<i,j> i 

The /ij's are random uncorrelated fields, distributed around zero, 

/li = 0, hihj = h^6ij, (3) 
We assume that the random fields are distributed according to a Gaussian distribution, 

P{h} = l[m,)^je-^^^^^\ (4) 

i 
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where A = {h\^)^. We are interested in calculating the average spin-spin correlations 
and susceptibility of a large, but finite, system, over a large number of realizations of the 
random field. Consider first 

rij = {aiaj)-{ai){aj). (5) 
The susceptibihty is related to the spin-spin correlations by 

^=|Er«- (6) 

The average susceptibihty, 

x = /3j]r^, (7) 
j 

is obtained by averaging over a large enough number of realizations. The true average of 
depends, of coarse, only on the radius separating i and j. Namely, translational invariance 
is restored by averaging. Obviously, since we consider a large system, calculating the above 
quantities directly involves the impossible task of performing a trace over a large number of 
spins numerically. We, thus, turn to real space renormalization. By choosing first a specific 
renormalization scheme (CS, MK, etc.), the rescaling factor, b, is set. We then choose the 
linear size of our system, L, to be b to some integer power, n >2, depending on how large we 
want it to be. The renormalization transformation is then used locally, by performing n — 1 
repeated iterations, to fully reduce the size of the system to b x b x b. The thermal average 
of each of the 8 remaining spins can now be calculated exactly by performing the trace using 
the Hamiltonian of the reduced system. In fact, averaging over large enough number of 
realizations, we may expect (cTj) to be translational invariant and, therefore, to be equal to 
the average magnetization per spin. In practice, in order to improve our averaging, we will 
calculate it as follows, 

_ 1 ^ 1 8 

M^j^j:(-^)n. = -sT.(-^)ns- (8) 
1=1 1=1 

Although, in principle, calculating thermal averages for the renormalized 8-spin system is 
a reasonable task, we are now faced with a different problem. The problem is that the 
reduced system only carries information about the 8 spins that survived the renormalization 
procedure. While this makes no difference when calculating quantities containing thermal 
averages of a single spin (such as the average magnetization above), it makes it impossible 
to calculate directly quantities that contain thermal averages of more than one spin. Such 
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are the average spin-spin correlations, Tij, at distances other than L/2, L/\/2 and ^/3L/2 
and such is the average susceptibihty, which, according to Eq. ((7j), requires the sum of Tij 
over all distances available in the original system. This is simply because spins that are 
nn's in the reduced system are, in fact, L/2 lattice constants apart in the original system 
(see Fig. H}. Indeed, as was done by Dayan at al [11], one may calculate the susceptibility 
by applying a small external field, H, to the original system and then use the derivative 
of M with respect to H. This method, though, is quite problematic, because it concerns a 
numerical derivative aX H = 0. As was discussed by Dayan at al, H must be small enough, 
so that M{H) is linear in H. This is difficult to achieve, since, below the transition, the size 
of the region where that linearity exists, shrinks to zero as the size of the system tends to 
infinity. If, on the other hand, the field is too small, one may encounter numerical problems 
from round-off errors. 

Our solution to the problem is obtained by using the identity: 



^ 1 



This identity was used in the past 



12, but since its proof is very short and simple. 



we will derive it here again for the sake of completeness of the presentation. Start with the 
right hand side of Eq. Q . 

1 ^-r^ 1 r, . d 



P{h}Dh. (10) 



1 fd{cx,) 



(3 J dhj 

This completes the proof since it is easy to see, using the random field Hamiltonian |21 that 

1 d{a,) 



(3 dhj 



{aiaj) - {ai){aj). (11) 



Connected spin-spin correlations and susceptibility. To obtain r( r ), the follow- 
ing procedure is used. The thermal average of a-i is calculated in a given realization for i that 
is one of the surviving spins, it is, then, multiplied by the value of hj in that realization, 
where j is a site separated by a vector from i on the original lattice. The product {<7i)hj 
is then averaged over many realizations. Since the true average should depend only on the 
vector connecting the sites and that up to a symmetry of the lattice, the statistics can be 
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FIG. 1: Starting with a large system of linear size L = 2"' (here demonstrated with n = 3), the 
system is fully reduced, using some RSRG transformation, to its minimal linear size of L = 2. 
The black numbered dots and the thickened lines connecting them, indicate, respectively, the 8 
remaining sites (only 6 of them are shown) and the remaining bonds connecting them, of the 
renormalized system. The circles and the doted lines reflect the boundary conditions imposed on 
the system. Sites that are nn's in the reduced system (such as 1 and 2) are, in fact, L/2 lattice 
constants apart in the original system. From all distances available in the original system (such 
as 3 <-> a and 3 <-> b), only L/2 for nn's, L/\/2 for nnn 's (like 1^4) and V3L/2 on the mam 
diagonals (like 3 6) are available in the reduced system. 

considerably improved by averaging r(^) as follows. 



where ik, runs over all equivalent sites around i, that are at a distance r from it. The number 
of these sites is rij.. In our numerical study, presented in Sec. IIVI we have limited ourselves 




(12) 
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to correlations in the directions of the principal axes of the lattice, so that rir — 6 (except, 
of coarse, for self correlations, where r — and rir — 1). For the average susceptibihty, the 
statistics is improved by writing 



^=^EwE^- (13) 
1=1 j 

In the following we will show how to evaluate "disconnected" correlation functions. Since, 
in practice, this involves much heavier computations, we just describe how it should be done, 
and postpone actual numerical application to future publications. 

The 'sites translation' method for the "disconnected" spin-spin correlations 
and the total average energy. The main point in the evaluation of "disconnected" 
correlations is that, actually, the method of integrating out many degrees of freedom and 
remaining with a small number of spins can yield not only ctj, where i is a surviving spin, 
but indeed all the local magnetizations. Namely, the method enables to calculate o"j for all 
i in the original lattice. This can be done by noting that, for a given realization, we can 
translate the 8 surviving spins. (This can be done, equivalently, by keeping the surviving 
spins and translating the field configuration.) Therefore, we choose a realization, evaluate 
the 8 ((Tj)'s corresponding to the surviving spins, then translate the field configuration by 
one lattice spacing, thus obtaining the local magnetization of the 8 sites translated from the 
original set by one lattice spacing in the opposite direction. This is repeated until all the 
local magnetizations are obtained. As is easily seen, this involves order of N repetitions of 
the original procedure. It is clear now how, by obtaining all the (cTj) 's, we can obtain 

1 



Sij = ^ij + = -^{'^i)hj + (14) 



This is generally time consuming but not so bad if we are interested, say, in the energy given 

by 

E=-J^J^-^1^, (15) 



<i,j> i 

for which our method yields the following expression 



<i,j> ^ i 

The calculation involves taking the ensemble average of the magnetization multiplied by 
the field at the same point (which is just Fjj) and the product where i and j ^-re 



nearest neighbors. To obtain the same degree of accuracy in Sij as in Tij, we need for each 
field configuration, to perform seven renormahzation procedures instead of the one needed 
for calculating Tij. The reason is that we need the original set of surviving spins and the 
sets obtained from it by the six unit translations in all directions. For other " disconnected" 
correlation, the time factor needed to attain the accuracy of the "connected" correlation is 

Hr + 1. 

Specific heat for systems of random bonds and fields. Next, we show how to 
obtain the average specific heat for an Ising Hamiltonian with both, random fields and 
random bonds. The evaluation of the average specific heat, C = kT"^ [O^"^) — (^)^j) is also 
made possible by our method. It is given by 



c " 



^{ij){ki) + 2 ^ hiJki^{i){ki) + ^ hihjVij j , (17) 

<k,l> <k,l> j / 



where for A and B, two sets of indices, we define Tab = {o'ao'b) — {o'a) (ctb) and = YlieA ^i- 
The parenthesis in the subscripts of V are used to describe the sets of spins. Operating on 
the three terms of Eq. ()17|) (from the less complicated on the right to the more complicated 
on the left), basically using integration by parts, we obtain first. 



hihjVij = hi (^hjVij + ^(^i)^ - ^hi{ai) 



d 

{hj{ai)) - {ui) 



dhj 



(18a, 



The next two terms are obtained similarly by following the same steps. The only difference 
is that, since we assume that the random bonds are distributed around some mean value, 
/ij, and uncorrelated with a standard deviation Aj, as a preliminary step, Jij and Jki are 
replaced with /ij + 5 Jij and /xj + 5Jki respectively. These are collected back by the end of 
the calculation. The results are similar to Eq. p8a|) only with the appropriate indices. 



hiJki^(i)(ki) = ^ 1^^^^ - 1^ hi{ai) (18b) 



and 



JijJki^(ij){ki) - ^-^^ - 1 - 5ij^ki^ Jij^OiOj). (18c) 

Now, Eq. ()18c|l is not yet in its final form, since it still contains the term (cjCTj). We thus 
first need to use Eq. ()14|) in order to fix that, and then Substitute Eqs. fjl8|) back into Eq. 
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It may seem plausible that techniques of the nature described above can be used also for 
correlations of quantities coupled by position independent coupling constants (that may be 
even zero). In principle, this is true because we can always add random couplings with a 
Gaussian distribution, do the calculation and, in the end, take the variance of those couplings 
to zero, practically, it is unclear whether such a procedure is more effective than taking a 
numerical derivative, since both imply taking the limit where certain couplings tend to zero. 
(In the random field problem, with no random bonds, we can either add random bonds with 
variance that will eventually tend to zero, or take a numerical derivative with respect to T.) 

IV. NUMERICAL RESULTS 

To demonstrate the usefulness of our method, we present in this section results for the 
less time consuming quantities, we present calculation of the average susceptibility and 
derive from it the critical exponents, comparing the results with those of Ref. jLll|, which 
uses the same renormalization scheme. We also present, for the first time we believe, a 
calculation of an V dependent correlation. We evaluate r(^) as a function of temperature 
and distance, for T^'s on the principal axes of the lattice. All figures are presented with 
error bars, although, in some of them, the error bars are too small to be noticed. The errors 
are the standard deviation calculated from the data. 

Using the CS approximation, we choose our input parameters, for calculating the average 
susceptibility, equal to those used by Dayan at al who also used the CS scheme. We 
also partially follow their line of analysis for extracting the critical exponents rj, 7 and u. 
This serves as a baseline, with which part of our results may be compered. We, thus, set 
J = 1 and h = 1 in the Hamiltonian Q and use as temperature. We have calculated the 
susceptibility as a function of temperature for systems of linear sizes L = 2,4, 8, 16, 32, 64 
and 128, averaging, each, over 10000 realizations of the random field, except for the largest 
system, for which we had to be satisfied with only 1000 realizations. This is shown in Fig. 
|21^a) in a Log-plot. Like Dayan at al, we also have included the mirror image of each random 
field realization, and that is in order to preserve, in our finite systems, the basic reflection 
symmetry of the infinite system. As seen in the figure, there is an upwards displacement 
of the average susceptibility for the largest, L = 128, system. This is probably because the 
relatively small number of = 1000 realizations, for that system, is too little to statistically 
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rely on. It may also be that, since the renormalization procedure in an approximation, the 
large number of renormalization steps generates an error which is too large. It is, therefore, 
excluded in the following analysis concerning the susceptibility. From Fig. El^b), we extract 
the value of 1], using the finite size scaling behavior of the susceptibility, 

X(T,(L)) ~ L^/^ (19) 

together with the scaling relation 

7 = (2-77)z/. (20) 

We thus obtain t] = 0.53 ± 0.003, where the error is the statistical one and errors generated 
by the renormalization approximate procedure are, therefore, not taken into account. 

To estimate the critical temperature, T^, of the infinite system together with the critical 
exponent 7, The logarithm of the average susceptibility of the largest system {L = 64) is 
plotted verses the logarithm of T — Tc, for temperatures above Tc and for different values 
of Tc. Acceptable values of T^. are such that, by lowering the temperature towards T^, 
the graph enters a linear region until finite size effects become important. We find that 
3.71 < Tc < 3.95. The two extremes are presented in Figs. Efa) andEfc). For critical 
temperatures above Tc = 3.95 the linear region disappears, while for critical temperatures 
around 3.71, an opposite curvature begins to appear, as demonstrated by Fig. Ofc) for 
Tc = 3.71. We consider Tc = 3.8 [Fig. E^b)], for which the largest linear region is obtained, 
to be the more probable value for the critical temperature. Next, the value of 7 for a given 
Tc is estimated using the critical behavior of the susceptibility, 

X(T)^A|T-Tc|-^ (21) 

where A is some constant. Taking the logarithm of both sides of Eq. 1)2111 . for a given Tc, 
the resulted straight line is, then, fitted to the linear region by varying 7 and the constant 
A. As indicated by Fig. El the range of acceptable Tc's corresponds to a range of possible 
values for 7: 1.9 < 7 < 2.4, which is quite similar to that obtained by Dayan at al jll|. 
Using the scaling relation ()20p . the resulting values for the critical exponent, u, are roughly: 
1.3 < ly < 1.65. For Tc = 3.8, the corresponding values for 7 and u are 7 = 2.2 and u = 1.5. 



In a search for a more refined estimation of gamma, we turn to the absolute value of 
the derivative of x with respect to T. Here we use the numerical derivative of our data for 
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Log[x(Tc{m 




(c) 



FIG. 2: In (a), the average susceptibility, x, is shown as a function of temperature, T, measured 
in units of 1//3. It is presented in a Log plot. The external field, H, is zero while the standard 
deviation of the random field is /i = 1. As indicated on each figure, the different levels of grayscaling 
correspond to systems of different linear size, L. While for L's up to 64, the system is averaged over 
10000 realizations, the largest, L = 128, system is averaged only over 1000 realizations. This may 
be the cause for the upwards displacement of the average susceptibility for the largest, L = 128, 
system. In (b), excluding the largest system, the logarithm of the maximum of the susceptibility 
is plotted verses the logarithm of L. The value of rj is obtained from a linear fit, as, by Eq. ()19() . 
the slope a is 7/1/ = 2 — t/. 
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FIG. 3: The logarithm of the average susceptibihty of a large system (L = 64) is shown as a 
function of the logarithm of T — Tc, for temperatures above Tc and for three different values of 
Tc. The two extreme, but yet acceptable, values of Tc, as discussed in the text, are presented in 
(a) and (c). It appears that the largest linear region is obtained for Tc = 3.8, as shown in (b). In 
(c), one may already notice that linearity breaks by the appearance of an opposite curvature. The 
value of 7 for a given Tc is determined by taking the logarithm of both sides of Eq. (|21() . 
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the average susceptibility, although one should note that it is, probably, better calculated 
directly from 

^ (22) 



dT 

3 

following our method and using the ideas presented in Sec. IIIII Similar to the finite size 
scahng behavior of x, as given by Eq. (|T9|l . its derivative with respect to T is expected to 
behave as 



dT 



(23) 



Here Tc{L) is the temperature of the maximum of the derivative for a given L and is, 
therefore, different from that of the maximum of used in Eq. p9|) . We have calculated. 



then, —dx/dT as a function of temperature for systems of linear sizes L = 2, 4, 8, 16, 32, 64, 
as presented in Fig. IH^a). In Fig. ^b), we have used a linear fit for the logarithm of the 



maximum of —dx/dT plotted versus the logarithm of L, from which we obtain the value of 
(7 + l)/z/. That, together with 7/1/ obtained earlier from the finite size scaling of x [Fig- 
I2fb)], determine the values of 7 and u. we obtain 7 = 2.41 ±0.07 and u = 1.64 ±0.07. Note 
that together with our previously obtained value of 77 = 0.53 ± 0.003, the scaling relation 
is satisfied. 

As a consistency check, we use the set of Tc(L)'s of the derivative (as they are better 
defined than those of the susceptibility itself), together with our estimation of z/ = 1.64, to 
extract the critical temperature Tc, directly from the finite size scaling behavior. 



|T,(L) - Tel ~ L-'/". 



(24) 



In Fig. El we use a linear fit for the logarithm of Tc{L) — plotted versus the logarithm 
of L. We obtain Tc = 3.77 and that is by tuning Tc to reach a slope that fits the negative 
inverse value of u = 1.64. 

Taking all the results above into account, we arrive at our final estimation of, 77 = 
0.53 ± 0.003, 7 = 2.2 ± 0.3, z/ = 1.5 ± 0.15 and Tc = 3.8 ± 0.1. Note that our result for rj 
satisfies the inequality 2 — 77 < d/2, bv Schwartz and Soffer jJJi] and is in good agreement 
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25|. Also, our result for 7 is in good 
27 L while our result for u is in good 



with the results cited in Refs. 
agreement with the results of Refs. 
agreement with the results of Refs 
to obtain a data collapse for the susceptibility, scaled by a factor of l/L^"'', when plotted 



29|. In Fig. 10 we use these values 
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(b) 



FIG. 4: In (a), the average of the absolute value of the derivative of the susceptibility is shown 
as a function of temperature. The external field, H, is zero while the standard deviation of the 
random field is /i = 1. As indicated on each figure, the different levels of grayscaling correspond 
to systems of different linear size, L. In (b), the logarithm of the maximum of —dx/dT is plotted 
verses the logarithm of L. The value of (7 + obtained from the slope of the linear fit here, 
together with 7/1^ obtained from the linear fit presented in Fig. Hfb), determine the values for 7 
and u. 
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LoglTciD- Tc] 




LogiD 



FIG. 5: The logarithm of Tc{L) — Tc, for systems of linear size L = 8, 16, 32, 64, is plotted versus 
the logarithm of L. The slope, a, is —l/u and Tc is tuned to 3.77 to make the slope of the linear 
fit to fit the negative inverse value of = 1.64. 
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versus T — T^, scaled by a factor of L~^/'^. It is presented in a log-log plot. It should be 
noted, though, that, within the range of values we have obtained for the parameters above, 
this data collapse picture, is almost insensitive, so that it is impossible to prefer one set of 
parameters over the other. 

We conclude our study of the average susceptibility by presenting it for different values 
of the strength of the random field (Fig. Ej). 

We now turn to our evaluation of the average spin-spin correlation function, F, calculated 



according to Eq. (fT^ as discussed in Sec. IIIII In Fig. |H1 we present F(r) for two values of 
the strength of the random field and for two temperatures. For h = 1, we know that T = 4.2 
is above the transition [Fig. |Hfa)]. We present our full results although it is clear that for 



r > 10 the values of r(r) are dominated by noise and therefore meaningless. For h = 2 
at T = 3.85 [Fig. lEfc)], we see a similar picture. Again the function decays very fast and 
already below r = 10, its significance is questionable. It may be expected that increasing the 
number of realizations considerably, may improve the evaluation of the correlation where it 
is small. As the temperature is lowered ,for h = 1, to T = 3.85 that is at the transition 
region [Fig. El^b)], the behavior becomes very noisy and statistically meaningless, but still a 
trend can be discerned. A similar behavior is observed in [Fig. Efd)], for h = 2 and T = 1. 
We have chosen to present Figs. lEfb) and Efd) although, as far as T(r) is concern ed, they 
are not very informative. The reason for doing so is that it is known |lll . 0, 0, 0, that 
self averaging is destroyed below the transition and thus we expect the enhanced noisiness 
of those figures to indicate approach to the transition. The existence of a trend suggests. 



though, that we may be still above the transition. Fig. |H1 presents r(r) for /i = 1 and T = 3 
(below the transition). The noisiness is larger than in Figs. El^b) andlHl^d) and no trend as 
a function of r can be observed. The above observations, in addition to information about 
the zero temperature transition 

n 

he = 2.27 according to [25|| and 



hr = 1.956 according to [2J|, while it is he = 2.28 and 
29^ respectively) is consistent with the qualitative phase 
diagram presented in Fig. ^| The above may suggest an alternative method of identifying 
the critical temperature by the amount of noise in the data and by loosing the trend as 
function of r. This line of investigation is postponed, however, to future work. 
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FIG. 6: The logarithm of the susceptibihty, scaled by a factor of l/L^^'', is plotted versus the 
logarithm of T — Tc, scaled by a factor of L^'^/^ . The data collapse is shown here for r/ = 0.53 
and u = 1.5, which, by Eq. (|2fl|) . corresponds to 7 = 2.2. The critical temperature is taken to be 
Tc = 3.8. 
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FIG. 7: The average susceptibility, x, is shown as a function of temperature, T. The external 
field, H, is zero, while, as indicated on each figure, the different levels of grayscaling correspond to 
different values of the standard deviation of the random field. The system is of linear size, L = 64. 
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FIG. 8: The average spin-spin correlation function, T, is shown as a function of the distance, r, 
taken along the main axes of the lattice and measured in units of lattice constant. The external 
field, H, is zero while the standard deviation of the random field is /i = 1 in (a) and (b) and h = 2 
in (c) and (d). As indicated on each figure, the two different levels of grayscaling correspond to 
systems of different linear size, L, indicating the size independence of F for these temperatures. 
The points of the larger system, with L = 128, appear to be more scattered since it is averaged only 
over 3500 realizations, while the smaller, L = 64, system is averaged over 10000 realizations. The 
different figures correspond to different temperatures. Note the broadening of F as the temperature 
is reduced towards entering the ordered phase at about T = 3.8 for /i = 1 and at T > for /i = 2. 
Also note that for h = 2, even at a temperature as low as T = 1, the correlations are kept relatively 
short-ranged indicating the persistence of the disorder phase to lower temperatures. 
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FIG. 9: The average spin-spin correlation function, F, is again shown as a function of the distance, 
r, As in Fig. |H1 only for a point on the phase diagram, (T,h) = (3.0,1.0), located below the 
transition. Note the noisiness and lost of trend compared with points in the phase diagram located 
above the transition [Fig. |S]. 
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FIG. 10: Schematic Phase diagram for the random field Ising system, are shown for d = 3 
(the thick fines). The zero temperature fixed point, {0,hc/J), controls the whole of the critical 
line, hc{T) [or Tc{h)], while the zero fields thermal fixed point, (Tc/J,0), that of the pure Ising 
system, is unstable (It is Tc/J = 4.57 for the Casher-Schwartz renormalization scheme 3). The 
lower horizontal dashed line, represent the lowering of temperature from the high temperature 
and disordered phase, at h = 1 < he- The critical line is, thus, crossed and the ordered phase 
is penetrated. In our simulations [Fig. [H], this is expressed by the fiattening of F. The higher 
horizontal dashed line, represent the same, only at h = 2 ^ he- As indicated by Fig. |H1 the level of 
flattening of F at a given point {T / J, h/J) on the phase diagram, depends on the distance of that 
point from the critical line. The X's in the figure represent the points used in the simulations. 
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V. SUMMARY 



We have presented a method for calculating thermodynamic quantities, directly from a 
fully reduced renormalized random system. Our method works with any renormalization 
scheme, though it is essential that the random variables are distributed according to a Gaus- 
sian distribution. It relies on an exact mathematical transformation, so that the quality of 
the results obtained by using it, depends solely on the quality of the renormalization approx- 
imation and the number of realizations considered. As examples, we have developed explicit 
expressions for the average "connected" spin-spin correlations, the average susceptibility, 
the average spin-spin correlations ("disconnected"), the total average energy and the aver- 
age specific heat. We have demonstrated our method by calculating the susceptibihty and 
the "connected" correlation function for the 3D random field Ising system. Prom the results 
for the susceptibility, we have calculated the following critical exponents: 77 = 0.53 ± 0.003, 
7 = 2.2 ± 0.3 and 1/ = 1.5 ± 0.15 while the critical temperature obtained is Tc = 3.8 ± 0.1 
for the case where the variance of the field is h — 1. As for the average "connected" spin- 
spin correlation function, we have presented it as a function of the distance between spins. 
Starting at the high temperature and disordered phase, it shows a sharp decay. Lowering 
the temperature towards the critical line it decays over longer and longer distances until the 
behavior becomes very noisy and no trend can be detected. 
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